function samples = gmmSampleSigmaPoints(mix, varargins)

if nargin < 3
    a = 1;
    b = 2;
    k = 0;
end

for i = 1:mix.ncentres
    samples(:, :, i) = gaussSampleSigmaPoints(mix.centres(i, :), mix.covars(:, :, i));
end



function samples = gaussSampleSigmaPoints(mu, sigma, varargin)

if length(varargin) < 3
    a = 1;
    b = 2;
    k = 0;
else
    a = varargin(1);
    b = varargin(2);
    k = varargin(3);
end

n = length(mu);
lambda = a*a*(n+k) - n;

[V D] = eig((n+lambda)*sigma);
scaledSqrtCovar = diag(sqrt(diag(D)))*V;
% scaledSqrtCovar = V*diag(sqrt(diag(D)))*V';

for i = 1:n
    samples(i, :) = mu + scaledSqrtCovar(i, :);
end
for i = n+(1:n)
    samples(i, :) = mu - scaledSqrtCovar(i-n, :);
end
% samples(end, :) = mu;
